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ABSTRACT 

We present the Millennium-II Simulation (MS-II) , a very large TV-body simulation 
of dark matter evolution in the concordance ACDM cosmology. The MS-II assumes 
the same cosmological parameters and uses the same particle number and output data 
structure as the original Millennium Simulation (MS) , but was carried out in a periodic 
cube one-fifth the size (100/i _1 Mpc) with 5 times better spatial resolution (a Plum- 
met' equivalent softening of 1.0 hr x kpc) and with 125 times better mass resolution (a 
particle mass of 6.9 x 10 6 M & ). By comparing results at MS and MS-II resolution, 
we demonstrate excellent convergence in dark matter statistics such as the halo mass 
function, the subhalo abundance distribution, the mass dependence of halo formation 
times, the linear and nonlinear autocorrelations and power spectra, and halo assembly 
bias. Together, the two simulations provide precise results for such statistics over an 
unprecedented range of scales, from halos similar to those hosting Local Group dwarf 
spheroidal galaxies to halos corresponding to the richest galaxy clusters. The "Milky 
Way" halos of the Aquarius Project were selected from a lower resolution version of 
the MS-II and were then resimulated at much higher resolution. As a result, they are 
present in the MS-II along with thousands of other similar mass halos. A compari- 
son of their assembly histories in the MS-II and in resimulations of 1000 times better 
resolution shows detailed agreement over a factor of 100 in mass growth. We publicly 
release halo catalogs and assembly trees for the MS-II in the same format within the 
same archive as those already released for the MS. 
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1 INTRODUCTION 

In order to understand how galaxies form and evolve in their 
cosmological context, we must understand the properties 
of dark matter halos over a wide range of physical scales 
and across virtually all of cosmic history. Numerical simula- 
tions provide one of the best methods for approaching this 
problem and have proven invaluable for studying the growth 
of cosmological structure and, in particular, of dark matter 
halos. Increasing computational power and improved algo- 
rithms have led to a steady and rapid increase in the abil- 
ity of A-body simulations to resolve the detailed internal 
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structure of dark matter halos over substantial cosmological 
volumes. 

Perhaps the most widely-used A-body simulation of 
cosmological structure formation to date has been the Mil- 
lennium Simulation (Springel et al. 2005 hereafter MS), 



which followed more than ten billion particles within a sim- 
ulation volume of (500 h -1 Mpc) 3 . This provided sufficient 
mass resolution to see the formation of halos hosting 0.1 L* 
galaxies and sufficient volume to obtain good statistical sam- 
ples of rare objects such as massive cluster halos and lumi- 
nous quasars. It also enabled the implementation of physical 
models for the formation and evolution of galaxy/AGN pop- 
ulations throughout a large and representative cosmological 
volume dCroton et al. 120061 IBower et al.||2006|). Since 2005, 



when the first results from the MS were published, most new 
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very large cosmological simulations have focused on larger 
volumes 1 ! (Lbox > 1000 hT x Mpc) in order to study topics 
such as the statistical detection of baryon acoustic oscilla- 
tions or weak lensing shear, or to build mock catalogs for 
the next generation of galaxy surveys ( Fosalba et al.|[2008 



Kim et al.|2008| |Teyssier et al.|2009| ). Moving to larger vol- 
ume simulations reduces computational cost at fixed par- 
ticle number both because resolved gravitational perturba- 
tions remain linear until later times and because the num- 
ber of simulation particles in a typical nonlinear structure 
is smaller. 

The opposite regime - smaller volumes with higher mass 
resolution - is much more computationally demanding but 
is also of great interest, especially for questions of galaxy 
formation, where the relevant mass scales are substantially 
smaller than for large-scale clustering. Understanding the 
formation and evolution of low-mass galaxies requires ade- 
quate resolution of the dark matter halos that host them, 
and this, in turn, requires much smaller particle masses 
than are currently feasible for Gigaparsec-scale simulations. 
These objects are important for galaxy formation as a whole 
because the first galaxies to form, which are low mass, pre- 
pare the initial conditions from which more massive sys- 
tems later form. Another topic of particular interest that 
can be addressed by high-resolution simulations is the evo- 
lution of substructure within dark matter halos. Such simu- 
lations show that subhalos can lose considerable mass after 
being accreted onto a larger halo - sometimes well over 99% 
without being completely disrupted ( Hayashi et al.| |2003 



Gao et al. 2004b Kravtsov et al. 20041. This means that 



as the resolution of a simulation is increased, so too is the 
typical time between accretion of a subhalo onto a larger 
system and its eventual tidal disruption. 

Moving to substantially higher resolution in a large- 
volume simulation is fraught with computational challenges, 
however. Increasingly small simulation time-steps are re- 
quired to accurately follow particle orbits in the dense cen- 



ters of dark matter halos ( Power et al. 2003 1, where the char- 
acteristic time-scale t s 



oc 1/y/G p is significantly shorter 
than on large scales. While only a small fraction of simu- 
lation particles reside in such dense regions, these particles 
are the limiting factor in how quickly the simulation can be 
evolved forward. The maximum resolved density contrasts at 
z < 1 can be one thousand times higher than those at z ~ 6; 
as a result, almost all of the computational time needed for 
such a simulation is spent at low redshift. Furthermore, the 
strong clustering of matter within a few very massive clumps 
can create serious problems with respect to parallelization: 
it is much more difficult to split such a particle distribution 
into optimal computational domains than is the case if the 
matter distribution is more homogeneous. 

In spite of these challenges, it is essential to have simu- 
lations that probe the structure of galaxy-scale dark matter 
halos with high mass resolution and over a large enough 
region to include a sizable and representative sample of 
objects. In this paper, we present such a calculation, the 
Millennium-II Simulation (hereafter MS-II). Section [2] gives 



details of the parameters which define this simulation and 
describes some of the post-processing we have carried out 
on its output, in particular, substructure-finding and merger 
tree-building. We present results on the evolution of the dark 
matter power spectrum and the two-point correlation func- 
tion in Section[3] In Section[4] we investigate the dark matter 
halo mass function and the clustering bias of dark matter 
halos. Section [5] focuses on halo formation times, including 
the dependence of clustering on formation time (so-called 
"assembly bias"; |Gao et al.|2005 I. A discussion of the rela- 
tion between the MS-II and the Aquarius Project ( |SpringeI] 
[et~aLl[2008| , as well as a comparison of the assembly histo- 
ries of the halos common to the two projects, is presented in 
Section [6] We summarize our results in Section [7] Through- 
out this paper, all logarithms without specified bases are 
natural logarithms. 



2 THE MILLENNIUM-II SIMULATION 
2.1 Simulation details 

The Millennium-II Simulation follows 2160' 3 particles within 
a cubic simulation box of side length Lbox = 100 Mpc. 
This is five times smaller than Lbox for the Millennium Sim- 
ulation. The volume sampled by the MS-II is thus 125 times 
smaller than in the MS but the mass resolution is corre- 
spondingly 125 times better: each simulation particle has 
mass 6.885 x 10 6 ft" 1 M Q . With this mass resolution, halos 
similar to those hosting Local Group dwarf spheroidals are 
resolved at our 20 particle mass limit, while halos of Milky 
Way-mass galaxies have hundreds of thousands of particles 
and halos of rich clusters have over fifty million particles. 
The Plummer-equivalent force softening] adopted for the 
MS-II was 1 h~ kpc and was kept constant in comoving 
units; this value corresponds to 0.06% (10%) of the virial 
radius for the largest (smallest) halos at redshift zero. 

The ACDM cosmology used for the MS-II is identical 
to that of the MS and the Aquarius simulations: 

fitot =1-0, fi m =0.25, n b = 0.045, n A = 0.75, 

h = 0.73, cr 8 = 0.9, n s = 1, (1) 

where h is the Hubble constant at redshift zero in units of 
100 km s _1 Mpc -1 , erg is the rms amplitude of linear mass 
fluctuations in 8 ft -1 Mpc spheres at z = 0, and n s is the 
spectral index of the primordial power spectrum. Retain- 
ing the cosmological parameters of the MS allows us to test 
for convergence by comparing results in the regime where 
objects are well-resolved in both simulations as well as to 
extend the range of structures probed by combining, when 
appropriate, results from the two simulations. In particu- 
lar, this helps us understand the effects of resolution at low 
particle number. 

The initial conditions for the simulation were created at 
redshift z — 127, identical to the starting redshift of the MS, 
using a "glass" initial particle load ( White|1996 1 ; the initial 
particle positions and velocities were then computed using 
the displacement field tabulated on a 4096 3 mesh and the 



1 Recent simulations with > 10 10 particles within smaller vol- 
umes (Lbox ~ 100 h _1 Mpc) have been used primarily for study- 
ing cosmic reionization at redshifts > 6 (e.g., Iliev et al.|2006 I. 



2 see eqn. 4 oflSpringell ||2005| and corresponding text for details 
of how force softening is implemented in GADGET. 
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Name 


Lbox 


Np 


e 




m p 


M min 


M max 


fgroup 




[h- 1 Mpc] 




[h^ 1 kpc] 


[h- 


-!M ] 


[h^ 1 Mq] 


[h- 1 Mq] 




Millennium-II 


100 


10,077,696,000 


1.0 


6.8 


9 x 10 6 


1.38 x 10 8 


8.22 x 10 14 


0.601 


Millennium 


500 


10,077,696,000 


5.0 


8.61 x 10 s 


1.72 x 10 10 


3.77 x 10 15 


0.496 


mini-MS-II 


100 


80,621,568 


5.0 


8.61 x 10 8 


1.72 x 10 10 


8.29 x 10 14 


0.502 



Table 1. Some basic properties of the new Millcnnium-II Simulation are compared to those of the MS and to the lower resolution 
version of MS-II (mini-MS-II). -Lbox is the side length of the simulation box, N p is the total number of simulation particles used, and e is 
the Plummer-equivalent force softening of the simulation, in comoving units. m p gives the mass of each simulation particle, M m j n gives 
the mass of the smallest FOF halos (corresponding to our choice of storing all halos with N p > 20), and Af max gives the maximum FOF 
halo mass found in the simulation. / gr0 up is the fraction of all simulation particles in FOF groups of 20 or more particles at z = 0. 



Zeldovich approximation. The transfer function used for cal- 
culating the input linear power spectrum was computed with 
the Boltzmann code CMBFAST ( jSeljak fc Zaldarriaga|1996| l. 

The amplitudes and phases of the initial linear fluctua- 
tion modes in the MS-II are identical to those in the simu- 
lation from which the Aquarius Project halos were chosen. 
Specifically, all modes with a wave-vector whose maximum 
component is less than 13.57 ft Mpc -1 have amplitudes and 
phases that match those of the Aquarius simulations; all 
other modes were set at random to have the same under- 
lying power spectrum. The Aquarius halos are thus present 
in the MS-II. A discussion of the relationship between the 
Aquarius Project and the MS-II is presented in Section [6] 

The MS-II was run with GADGET-3, an updated version 
of the GADGET code ( |Springel et al.||2001b| |Springel|[2005j ). 
GADGET-3 is a TreePM code: long-range force calculations 
are performed with a particle-mesh algorithm while short- 
range forces are calculated via a hierarchical tree. While the 
original MS was performed with a memory-optimized ver- 
sion of GADGET-2, the extremely high level of clustering that 
occurs in the MS-II results in somewhat different compu- 
tational requirements; in particular, a more flexible domain 
decomposition is necessary. GADGET-3 was developed specif- 
ically for this situation. 

The MS-II was performed on the IBM Power-6 com- 
puter at the Max-Planck Computing Center in Garching, 
Germany, using 2048 cores and approximately 8 TB of mem- 
ory. A Fast Fourier Transform with 4096 3 cells was used for 
the PM calculation. Particles were allowed to have individ- 
ual, adaptive time-steps. The evolution of the simulation re- 
quired approximately 1.4 million CPU hours and 2.77 x 10 13 
force calculations for the 22,142 simulation time-steps. The 
mid-point of the simulation in terms of computational time 
was z = 0.88; by contrast, evolving the simulation from 
z = 127 to z = 6 took only 10% of the total CPU time. 

Outputs were saved at sixty-eight epochs: sixty-five 
snapshots spaced according to 

log 10 (l + ZJ v)= N( " 2 Q 35) (0 <iV<64) (2) 

and three high redshift outputs at z = 40, 80, and 127 . The 
spacing scheme in Equation Q is identical to that used in 
the MS. We have extended the range of regularly-spaced 
outputs to z m 31.3, however, because the increased mass 
resolution of the MS-II results in earlier-forming first struc- 
tures. 

For comparison purposes, we have also performed a ver- 
sion of the MS-II with identical initial conditions and in the 
same volume with the same outputs as the main run but 



at the same mass and force resolution as the original MS 
(so N p = 432 3 ). This ' 'mini-MS-II" simulation allows us to 
test how numerical resolution affects our results. Some basic 
details of all three simulations are listed in Table [T] 



2.2 Halos and Subhalos 

Dark matter halos were identified on-the-fly during the sim- 
ulation for each snapshot using the friends-of-friends (FOF) 
algorithm ( Davis et al. 1985 1 with a linking length of b = 0.2; 



all groups with at least 20 particles were retained for later 
analysis. This process resulted in 1.17 x 10 7 FOF groups at 
2 = 0, slightly fewer than the peak value of 1.53 x 10 7 at 
z = 3.06. Just over 60% of the particles in the full simulation 
belong to a FOF group at z — 0. A catalog with quantities 
of interest for each FOF halo (e.g., position, velocity, num- 
ber of particles), as well as a list of the particles composing 
each halo, was saved at each snapshot. 

The largest FOF group at z — 0, a cluster-mass dark 
matter halo, has over 119 million particles. Figure [T] shows 
images of the dark matter distribution in the MS-II on a 
number of different physical scales, all centered on this haltj^J 
The large panel in the upper left shows a 15 ft -1 Mpc-thick 
slice through the full simulation volume (100 h~ Mpc on 
each side). The well-known cosmic web of filaments and 
voids can be seen clearly. Starting in upper-right and mov- 
ing clockwise, the other five panels zoom successively closer 
into the halo. The bottom-right panel is 5 h~ l Mpc on a side, 
approximately the diameter of the halo. As has been long 



known (e.g., Moore et al. 


1998 


Tormen et al.|1998 Ghigna 


et al.||1998 Klypin et al. 


1999a 


b Moore et al.|1999l, FOF 



halos in ACDM simulations are not monolithic objects but 
rather are teeming with substructure; this substructure is 
clearly evident even at l/10 th the radius of the halo (lower- 
left panel). 

During post-processing, every FOF halo was searched 
for bound dark matter substructure using the SUBFIND al- 
gorithm ( Springel et al.|2001a |. SUBFIND identifies substruc- 



tures within a FOF halo by searching for overdense regions 
using a local SPH density estimate, identifying substruc- 
ture candidates as regions bounded by an isodensity surface 
that traverses a saddle point of the density field, and testing 
that these potential substructures are physically bound with 
an iterative unbinding procedure. All self-bound structures 

3 Images of individual panels and additional information related 
to the MS-II are available at 



http: / /www.mpa-garching.mpg.de/galform/millennium-II 
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Figure 1. A sequential zoom through the Millennium-II Simulation. The large image (upper left) is a 15 h~ 1 Mpc thick slice through 
the full 100 h^ 1 Mpc simulation box at redshift zero, centere d on the most massive halo in the simulation. This FOF halo has Mfof = 
8.2 X 10 14 h~ 1 Mq, similar to the mass of the Coma cluster |Colless & Dunnll996k, is composed of 119.5 million particles, and contains 
approximately 36,000 resolved subhalos spanning 6.7 decades in mass. Starting from the upper right and moving clockwise, subsequent 
panels zoom into the cluster region and show slices that are 40, 15, 5, 2, and 0.5 h' 1 Mpc on a side (with thicknesses of 10, 6, 5, 2, and 
0.5 h^ 1 Mpc). Even at 0.5 h~ 1 Mpc, which is approximately l/10 th the diameter of the halo, a rich variety of substructure is visible. 
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Figure 2. Time evolution of the largest FOF halo at z = in the Millennium-II Simulation. The halo is shown at three co-moving scales 
(from left to right: 100, 40, and 15 h~ 1 Mpc, with thickness 15, 10, and 6 Mpc) and at four different cosmological epochs (from top 
to bottom: 2=6.2, 2.07, 0.99, and 0). 
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with at least 20 particles were deemed to be physical subha- 
los and were stored in subhalo catalogs. Several properties 
of each subhalo were also tabulated and saved, including ve- 
locity dispersion, peak circular velocity Kn ax and the radius 
-Rmax at which Vmax is attained, half-mass radius, spin, po- 
sition, and velocity. The member particles of each subhalo 
were ranked according to binding energy and stored in that 
order, which facilitates tracking subhalos across simulation 
outputs. Note that with these procedures, we have two sep- 
arate but related sets of dark matter structures: FOF halos 
and subhalos. 

While each subhalo has a single well-defined mass as- 
signed to it - the sum of the masses of its constituent parti- 
cles - multiple mass definitions for FOF halos are common 
in the literature (see [White 2001 for a discussion of sub- 
tleties associated with assigning masses to halos) . The most 
straightforward definition is Mfof, the total mass of all the 
member particles. Another possibility is Ma, defined as the 
mass contained in a spherical region (centered on the parti- 
cle in the dominant subhalo with the minimum gravitational 
potential) with average density a factor A larger than the 
critical density of the universe. For each FOF halo, we cal- 
culated M A for A = 200, 200f2 m (z), and A v (z), where the 
last value is taken from the spherical top-hat collapse model 
(see, e.g., Bryan fc Norman|1998[ ). We refer to spherical over- 
density masses as M 20 o [A = 200], M 20 om [A = 200fi m (z)], 
or M v [A = A„ (z)] and to the corresponding virial radii as 
R200, fl^oom, or R v . At high redshifts, when the matter den- 
sity is nearly equal to the critical density, all three definitions 
give similar masses. At lower redshifts, M2oom > M v > M200 
for a given halo. 



2.3 Merger Trees 

Merger trees were constructed at the subhalo level by re- 
quiring subhalos to have at most one descendant. For many 
subhalos, this descendant can be found trivially (if it exists): 
all particles in a subhalo at snapshot S n may belong to a 
single subhalo at the subsequent snapshot S n +i, in which 
case this subhalo is clearly the descendant of the subhalo 
at the previous snapshot. There is also the possibility that 
particles belonging to one subhalo at S n may be distributed 
over more than one subhalo at SVi+r- We still require each 
subhalo to have at most one descendant for these cases, so 
a subhalo's unique descendant is identified as follows. First, 
the binding energy of each particle in the subhalo at S n 
is calculated and the particles are ranked by this binding 
energy. Each potential descendant subhalo - that is, each 
subhalo at S n +i containing at least one particle j from the 
subhalo at S n - is then given a score \ that is based on the 
binding energy rank 1Z of these particles: \ ~ ^2jT^-J 2 ^ 3 ■ 
The subhalo at S n +i with the largest value of \ is defined to 
be the descendant subhalo. This procedure weights the most 
bound regions of a subhalo most heavily when determining 
its descendant. Note that while descendants are unique, a 
given subhalo may have many progenitors. 

There is one slight complication to this process. Some- 
times a subhalo passing through the dense center of a larger 
system will not be identified by SUBFIND, simply because 
the density contrast is not high enough. To mitigate this 
problem, we also search for a descendant at snapshot S n +2- 



In the vast majority of cases, however, the descendant of a 
subhalo is found at S n +i- 

Once all unique descendants are found, the subhalos 
are linked across all snapshots to form merger trees. This is 
done by taking a subhalo at z = and linking all subhalos 
with descendant pointers to this halo, then repeating with 
all of those subhalos, and so on, until no more subhalos 
can be joined. This process results in links between most, 
though not all, of the subhalos in the simulation: subhalos 
that are never connected to any z = subhalo and that 
are never connected to any progenitor of any z — subhalo 
are not included in the trees. We save several pointers for 
each tree subhalo for later use. These include pointers to 
the dominant subhalo of the subhalo's FOF group, the next 
most massive subhalo in the FOF group, the progenitor that 
contains the largest fraction of the subhalo's particles, the 
subhalo's descendant, and the next most massive subhalo 
that shares the same descendant 

The merger trees for the MS-II contain approximately 
590 million subhalos in total (as compared to 760 million 
subhalos in the MS). While the overall data volume of the 
MS-II is similar to that of the MS (« 25 Terabytes, domi- 
nated by the raw particle data), the highly clustered nature 
of the MS-II means that the trees are markedly less homo- 
geneous. There are only half as many trees in total in the 
MS-II as in the MS, but the largest tree is much larger, with 
over 90 million subhalos (compared to 500 thousand for the 
largest tree in the MS). 



2.4 An example of subhalo tracking 

As an example of our merger trees and subhalo tracking, 
we consider the main progenitor histories of the most mas- 
sive halo at z = and of two of its subhalos. We use both 
the MS-II and mini-MS-II in order to highlight the effects 
of resolution and to probe the convergence of the subhalo 
identification and merger tree building algorithms. Objects 
can be matched between the simulations because the initial 
conditions are identical on all scales that overlap; any differ- 
ences are due to force and mass resolution and to differing 
discreteness effects. 

The main cluster halo is trivial to find in both simula- 
tions. At z = 0, the properties agree quite well between the 
two: the FOF masses are the same to within 1% (Mfof = 
8.22 x I0 15 h" 1 M e for MS-II versus 8.29 x 10 15 h" 1 M for 
mini-MS-II) and even the position of the halo, as deter- 
mined by the gravitational potential minimum, agrees to 
with 0.022 ft" 1 Mpc = 4e (for mini-MS-II). Figure [| shows 
the main progenitor of this halo on three co-moving scales 
(from left to right: 100, 40, and 15 Mpc) and at four 
redshifts (from top to bottom: 2 = 6.2, 2.07, 0.99, and 0.0). 
Using the merger trees, we track the main progenitor of the 
most massive cluster back in time for each run until there 
are no further progenitors. The top two curves in Figure [3] 
show the mass of the central subhalo of this main progenitor 
branch^ for the MS-II (black) and mini-MS-II (magenta). 



see also figure 5 in the Supplementary Information of |Springel| 
et"al] ( [2005) 

u Note that the while the FOF masses of the halos agree to within 
1% between the MS-II and mini-MS-II, the mass of the central 
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Figure 3. Subhalo assembly histories in the MS-II and mini-MS- 
II. The upper set of curves shows the mass in the main progenitor 
branch of the most massive 2 = FOF group for the MS-II (black) 
and mini-MS-II (magenta), while the lower set of curves shows 
main progenitor histories for subhalos that are accreted onto the 
FOF group of the cluster at z ~ 2 (solid lines) and z 3 (dashed 
lines). The horizontal dotted line shows the mass resolution of 
MS and mini-MS-II, while the vertical dotted and dot-dashed 
lines show the epochs at which the smaller halos joined the FOF 
group of the main cluster in the MS-II. 

The two are in excellent agreement from z — to z = 9, 
at which point the main branch from mini-MS-II falls be- 
low 100 particles and resolution effects become relevant (the 
main progenitor in the MS-II can be traced back all the way 
to z « 18). Clearly, the assembly of the main progenitor is 
very well converged between the two runs at z < 9. 

We also consider the evolution of two far less massive 
subhalos within this FOF group: subhalo A (lower solid lines 
in Figure |3| and subhalo B (dashed lines) are identified 
in the MS-II at z = 0. They are not massive enough to 
be identified at z = in mini-MS-II, as subhalo A has 57 
particles and subhalo B has 812 at the final snapshot (the 
resolution limit of mini-MS-II corresponds to 2500 particles 
from the MS-II and is marked by the horizontal dotted line). 
Nevertheless, Figure [3] shows that by tracking subhalos A 
and B backward (lower black curves), we find that both 
were much more massive in the past: each exceeded 5 x 
10 11 ft" 1 M© at one point in its history (92,431 particles for 
A and 73,718 for B) and each was the central subhalo of its 
own FOF group before falling into the main progenitor of 
the cluster (this accretion happened at z ~ 2 for subhalo A, 
marked by the vertical dotted line, and at z ~ 3 for subhalo 
B, marked by the vertical dot-dashed linej^] 

subhalo in the MS-II is slightly smaller because more distinct 
subhalos are identifiable. Many subhalos that are resolvable in 
the MS-II but not in mini-MS-II show up as extra mass in the 
central subhalo in mini-MS-II. 

6 A and B are two of four subhalos in the main FOF cluster 



Figure 4. The mean cumulative subhalo abundance per parent 
halo in the MS-II (solid curves) and mini-MS-II (dashed curves) 
in five parent halo mass bins. The curves for each simulation are 
plotted down to the 20 particle resolution limit. There is a deficit 
in subhalos at this limit in mini-MS-II relative to the MS-II due to 
resolution. At 3-5 times the minimum resolution limit, however, 
the two simulations agree very well, indicating that subhalos with 
more than 100 particles are reliably resolved. 

These maximum masses for A and B are easily resolv- 
able in mini-MS-II, so we can hope to find the equivalent 
subhalos there and compare their mass histories. It is in 
general unrealistic to expect subhalos to have identical po- 
sitions in runs of differing resolution: the gravitational force 
is softened at different scales and there is a difference in the 
'graininess' of the gravitational potential due to finite parti- 
cle mass, both of which can cause orbital phase offsets that 
accumulate over tim^] Nevertheless, we are able to locate 
A and B in mini-MS-II and we find them to be at almost 
exactly the same positions as in the MS-II: the absolute po- 
sitions for subhalo A (subhalo B) differ by only 0.015 (0.010) 
hT 1 Mpc at the times marked by the vertical lines, which is 
only three times the force softening of the low resolution 
run. 

Having located A and B in mini-MS-II, we then use the 
merger trees to track the subhalos both forward and back- 
ward in time and we compare to the results from the MS-II. 
These subhalo mass histories are shown in the lower ma- 
genta lines (solid for subhalo A, dashed for subhalo B) of 

at 2 = that (i) have N p (z = 0) < 1000 and (ii) have over 
1500 progenitor subhalos in their sub-tree. This selection picks out 
halos that were massive at one point in their history but are not 
at 2 = 0. We have also investigated the other two subhalos from 
this sample and find similar convergence in the subhalo tracking 
and mass identification between the two simulations; for clarity, 
the res ults are not plo tted in Figure [3] 

7 See |Springel etaT] j2008| l for a method to match subhalos in 
simulations based on the positions of the subhalo particles in the 
initial conditions. 
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Figure [3] We see that the subhalos in the two simulations 
have remarkably similar assembly histories, not only when 
they are the main FOF subhalo (to the right of the vertical 
lines) but also when they are non-dominant subhalos within 
a larger FOF group (to the left of the vertical lines). This 
regime, where the subhalos are subjected to strong tides that 
vary rapidly in time, can be extremely difficult to capture 
accurately in simulations of differing resolution. The excel- 
lent agreement between the lower magenta and black curves 
demonstrates that the subhalos have the same dynamical 
histories, are assigned the same masses, and are linked in 
the same way by the merger trees in the two runs. 

Figure [3] also illustrates resolution limitations. As tides 
strip mass from the subhalos, they are lost from the mini- 
MS-II catalog and are considered to have merged with the 
dominant subhalo, while they persist as independent subha- 
los to z = at the significantly enhanced resolution of the 
MS-II. The high resolution of the MS-II is required to study 
the fates of subhalos hosting low mass galaxies within larger 
structures: note that the maximum masses - approximately 
6 x 10 11 hT 1 Mq - of A and B are quite large, larger than the 
halo masses of likely Milky Way progenitors at the redshift of 
accretion into the massive halo. A and B are therefore likely 
to host galaxies of stellar mass comparable to that of the 
Milky Way. At mini-MS-II resolution (i.e. MS resolution) it 
is not possible to follow the dynamics of these subhalos past 
z w 1.5 for subhalo B or z ~ 0.5 for subhalo A, by which 
redshifts their masses have dropped below 10 11 h~ x Mq. The 
MS-II captures the later dynamical history of both subha- 
los, even though subhalo A (subhalo B) retains only 0.06% 
(1.1%) of its maximum mass at z = 0. Note that these fi- 
nal masses are considerably smaller than the likely stellar 
masses of the associated galaxies, so it remains unclear how 
realistic the late-time dynamical evolution actually is in MS- 
II. 

We can also consider the statistical agreement between 
the MS-II and mini-MS-II by comparing stacked subhalo 
abundances as a function of host halo mass. Figure [4] shows 
the mean number of subhalos per host halo in five host halo 
mass bins for the two simulations. Resolution effects reduce 
the number of subhalos at a given subhalo mass in mini- 
MS-II (dashed curves) relative to the MS-II (solid curves) for 
subhalos with few particles: at the minimum resolvable mass 
of 20 particles in mini-MS-II, the abundance of subhalos is 
reduced by approximately 30% relative to the MS-II. The 
results from the two simulations are in excellent agreement 
for more massive subhalos (Af su b > 10 11 Mq), showing 
that subhalos containing at least 50-100 particles are reliably 
resolved. 



3 STATISTICS OF THE DENSITY FIELD 
3.1 Power Spectrum 

At comoving position x and time t, the mass density field 
can be expressed as 



p(x,t) =p(t) [1 + J(x,i)]. 



(3) 



In the standard picture of structure formation in a cold dark 
matter universe, the initial density fluctuation field 5(x, 0) is 
taken to be a Gaussian random field. Its statistical proper- 
ties are therefore fully specified by its power spectrum P(k) 




0.8 



1 10 
k [hMpc- 1 ] 



1000 



Figure 5. The power spectrum A 2 (fc) measured from the MS-II 
at redshifts 0.0 (black curves), 0.99 (magenta), 2.07 (green), and 
6.20 (cyan), as well as the linear theory power spectrum at each 
redshift (gray curves). Power spectra from the MS (dotted curves) 
at the same redshifts are also shown for comparison. The dashed 
lines correspond to the shot noise limit for the MS-II (right) and 
the MS (left); the power spectra have not been corrected for shot 
noise. The bottom panel shows the ratio of the power spectra. 



or equivalently its dimensionless power spectrum A (k). 

k 3 



2tt 2 



P(k). 



(4) 



A 2 (fc) measures the power per logarithmic interval in 
wavenumber; A 2 (fc) ~ 1 therefore indicates that fluctua- 
tions in density on scales near wavenumber k are of order 
unity. 

The dark matter power spectrum from the MS-II is 
shown in Figure [5] We plot the results at four redshifts: 
z = (black curves), 0.99 (magenta), 2.07 (green), and 6.20 
(cyan). On large physical scales (small wavenumber k), the 
power spectrum follows the prediction from linear theory 
(light gray lines). As time progresses, larger and larger phys- 
ical scales become non-linear and the small-scale power ex- 
ceeds the linear theory prediction. Results from the MS are 
also included in Figure [5] for comparison (dotted curves). 
The agreement between the two simulations is very good, 
and the MS-II extends the measurement of A 2 (fe) by a fac- 
tor of 5 at large k. We have not performed a shot noise 
correction in this figure, as it not clear that it is appropriate 
to do so (see, e.g., |Baugh et al.||1995| |Sirko|2005| . The shot 
noise limit for each run is plotted as a dashed gray line. 

The large dynamic range and uniform mass resolution 
of the MS-II allows us to probe the dark matter power spec- 
trum on a wide range of scales, including scales where ex- 
isting fitting functions (Peacock & Dodds 1996; Smith et al. 
2003 1 are uncalibrated and untested. Figure |6| compares the 
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Figure 6. A comparison of the MS-II power spectra at four red- 
shifts to the halo model fit from Smith et al. l |2003| gray curves). 
The fit at z = is accurate for k < 7 h Mpc -1 but underestimates 
the power spectrum from the MS-II by 50% or more for k between 
10 and 100 h Mpc -1 . Redshifts one and two show similar results 
but the agreement at z = 6 is poor. 



power spectrum computed from the MS-II with the halo 
model predictions of Smith et al. (gray lines). At redshift 
2 and below, the Smith et al. model agrees with the cal- 
culated power spectrum to within 10% for k < 5 ft Mpc . 
At larger k, however, the model significantly underestimates 
the power, with errors exceeding a factor of 2. The Smith 
et al. model was not calibrated for this range, so it is not 
surprising that it fails to reproduce the simulation results. 
Nevertheless, this is a reminder that extrapolating fitting 
functions beyond their calibrated range can lead to serious 
errors. At z — 6, the Smith et al. model fails to fit the data 
over the range 1 < k < 10 h Mpc" 1 , where the MS-II and 
MS agree very well. 

At scales of k > 10 h Mpc - , baryonic physics plays an 
important role in determining the real dark matter power 
spectrum (e.g., Rudd et al.|2008| ). Although a full modeling 
of baryonic effects will be necessary to get accurate predic- 
tions for A 2 at these scales, understanding the underlying 
dark matter-only power spectrum still provides a critical 
baseline. 



3.2 Two Point Correlation Function 

The spatial two-point correlation function of the density 
field is given by 



Figure 7. Measurements of the two point correlation function 
§ as a function of comoving separation r from the MS-II. Wc 
show four redshifts: z = 0.0 (black curves), 0.99 (magenta), 2.07 
(green), and 6.20 (cyan) and we compare with §(r) from the MS 
(dotted curves) at the same redshifts. On large scales, the correla- 
tion functions from the two simulations agree quite well. On small 
scales (< 0.020 h -1 Mpc in physical units), the MS-II correlation 
function amplitude is larger, reflecting structures that are not re- 
solved in the MS. The bottom panel focuses on these differences 
by plotting the ratio of the correlation function from the MS-II 
to that from the MS. 



or equivalently, by a Fourier transform of the power spec- 
trum: 



<;< )- / A 2 (fc)^Mdlogfc 



(6) 



The correlation function is a useful measure of the spatial 
clustering of dark matter: it gives the excess probability of 
finding pairs of particles at a given separation relative to 
a Poisson distribution. Figure [7] shows £(r) at redshifts 0, 
0.99, 2.07, and 6.20, with results from the MS at the same 
redshifts also plotted for comparison. (Note that the scale 
on the horizontal axis of Figure [7] is in co-moving units.) 

The correlation function shows a prominent feature at 
r « 2ft" 1 Mpc (for z = 0). This is the well-known transi- 
tion between the 'one-halo' and 'two-halo' contributions: on 
smaller scales the correlation function is dominated by dark 
matter particle pairs within the same halo, while on larger 



scales it is dominated by pairs in separate halos (Peacock & 



Smith 2000; Scljak 2 000] |Ma fc Fry|2000| |Scoccimarro"et~aT 
2001j|Cooray fc Sheth |2002[). At no redshift is the correlation 



£(r) = (<5(x)«5(x + r)> 



(5) 



function even roughly approximated by a single power law. 
This is in stark contrast to observations of galaxy correla- 
tion functions ( jZehavi et al.|2 002 Hawkins ct al. 200 3| and 
the stellar mass autocorrelation function ( Li fe White|2009| 
at low redshift, which show a remarkably good power-law 



© 2009 RAS, MNRAS 000,[T|fl6] 



10 M. Boylan-Kolchin et al. 



behavior £ oc r~ 1,8 over 1CP 2 < r < 10 hT x Mpc (a gray line 
with this relation is also plotted in Figure [7] for comparison). 

From z = 2 to z = 0, the agreement between the MS-II 
and MS results is quite good from 0.03 to 2 ft -1 Mpc. The 
MS-II result lies approximately 10% below the MS £(r) on 
scales of 2-10 h~ x Mpc; this is the range where the contribu- 
tions from pairs in separate halos become important and is 
likely an indicator of cosmic variance in the two-halo term. 
On small scales, the MS-II correlation function lies above 
that of the MS. This is a result of the higher force and mass 
resolution of the MS-II: low-mass halos that were not re- 
solvable in the MS also boost the clustering on small scales. 
The MS correlation function is noticeably lower in amplitude 
than £(r) from the MS-II at z = 6 for r < 0.2 h' 1 Mpc. This 
coincides with the mean interparticle spacing, 0.231 ft -1 Mpc 
comoving for the MS, and is therefore most likely due to dis- 
creteness in the glass-like particle load used for the initial 
conditions. 



4 DARK MATTER HALOS 

Understanding how dark matter overdensities grow and ul- 
timately virialize into highly non-linear structures is an ex- 
tremely difficult problem from a theoretical perspective. No 
rigorous analytic techniques are available for use in both the 
linear and non-linear regimes. The most successful model for 
dark matte r halo formation ([Press fc Schechter]|1974 Bond 



10.0 



et al.|1991| |Lacey fc Cole|1993| |Sheth et al.|2001| see [Zent 
ncr 1 120071 for a recent review of extended Press-Schechter 
theory) relates the abundance of halos at mass M and red- 
shift 2 to the initial power spectrum of density fluctuations 
and to the well-understood regime of linear growth. It con- 
siders the linear overdensity field smoothed using a spherical 
top-hat filter (in real space) containing mass M and extrap- 
olated using linear theory to redshift z. The variance of this 
smoothed field is: 

(j 2 {M,z) = d 2 {z) J A? in (fc) W 2 (k;M)dlogk , (7) 

where d(z) is the linear growth factor at redshift z with 
normalization d(z = 0) — l, Af in (fc) is the linear power spec- 
trum extrapolated to z = 0, and W(k; M) is the Fourier 
transform of a real-space spherical top-hat filter. Each dark 
matter particle is assigned to a halo of mass M at redshift 
z, where M is taken to be the largest filter scale for which 
the smoothed linear overdensity at the particle's position 
(extrapolated to redshift z) exceeds a threshold value <5 C [^] 
A characteristic halo mass can then be defined at each 
redshift via a(M±, z) = 8 C . In this model, the halo multiplic- 
ity function /(<r), which is related to the comoving number 
density of halos via 



M 



dn(M, z) 

d log cr _1 



(8) 



takes on a universal form (for perspectives on universality of 
the halo mass function, see Jenkins et al. 2001; Reed et al. 



8 <5 C = 1.68647 at all redshifts for an Einstein-de Sitter universe. 
For our ACDM cosmology, S c varies slightly, from this standard 
value at high redshift to 1.6737 at redshift zero due to a weak 
dependence on Q m (z). 
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Figure 8. The halo multiplicity function f(cr) as a function of 
cr(M, z)~ 1 at four redshifts from the MS-II: z = (black squares), 
0.99 (magenta), 2.07 (green), and 6.2 (cyan). We compute the 
multiplicity function both with the Warren et al. N p correction 
(lower data points) and without the correction (upper points, 
offset by 1 dex for clarity). Also overplotted is the [Warren et ah] 
p006\ fit to the halo multiplicity function. The MS-II multiplicity 
function shows universal behavior when scaling with respect to 
redshift, with deviations at the 10% level. Halo masses here are 
defined to be MpQF- 



2007 



2008) 



|Lukic et al.|[200"7l |Cohn fc Whltel[2008l [Tinker et al.| 



4.1 Mass Function 

In Figure[8]we plot the halo multiplicity function f(a) from 
the MS-II at four redshifts: z = 0.0, 0.99, 2.07, and 6.20. 
Here we define halo mass as Mpof; see |Hu fc Kravtsov] 
( |2003[ ) and |Tinker et al.| ( |2008| for discussions of halo 
multiplicity functions computed using spherical overdensity 
masses. We include Poisson error bars for all bins contain- 
ing fewer than 400 halos, as Poisson errors dominate sample 
variance at all masses (e.g .,|Hu fc Kravtsov||2003[ ). Results 
are plotted both using the Warren et al.| ( |2006[ ) correctiorp] 
for sampling bias in N p (lower set of data points) and with- 
out this correction (upper data points, offset upward by 1 
dex for clarity). We exclude halos that have N p < 20 when 
using the corrected N p . The simulation has fixed mass res- 
olution but a(M,z) evolves significantly with time, so by 
comparing multiplicity functions at several redshifts we can 
probe a large range in a. 

The multiplicity function within the MS-II does seem to 
have a universal form: where the data overlap, the agreement 
in f(a) is quite good. It thus appears possible to compute 
the multiplicity function at any redshift simply by combin- 
ing the linear growth factor d(z) with the rms amplitude 



This correction is N p — > iV correc t e< i = N p (1 — N p °' 6 ) 
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Figure 9. FOF mass function for the MS-II (solid blue squares) 
and for the MS (open red squares) compared at redshift 6.2 and 
0. The redshift zero mass functions are in excellent agreement 
over the entire range where the two simulations overlap. At red- 
shift 6.2, the MS-II points lie systematically above those from the 
MS. The shaded gray region shows the range of mass functions 
obtained from subdividing the MS into 125 cubes with volume 
equal to the MS-II and computing a mass function for each sub- 
volume. The MS-II points are well within the scatter, indicating 
that the difference is likely due to the small volume of the MS-II. 



of fluctuations as a function of mass at redshift zero. This 
agreement is at least as good for the uncorrected points, ex- 
cluding bins containing halos with fewer than 100 particles. 
(We note, however, that Warren et al. used a minimum of 
400 particles per halo in deriving their fitting parameters; in 
this regime, both the corrected and uncorrected points seem 
to exhibit 'universality.') The multiplicity function does not 
agree precisely with the Warren et al. fit (gray line) in ei- 
ther case; however, the volume of the MS-II is not sufficiently 
large to obtain statistically precise results in the high a" 1 
regime due to cosmic variance. 

Figure [9] compares the FOF mass function at redshifts 
and 6.2 determined from the MS-II (solid blue squares) 
with the MS mass function (open red squares). Poisson er- 
ror bars are included for all bins with fewer than 400 halos 
and the data points do not include the Warren et al. correc- 
tion for the sampling bias in N p . At z = 0, the agreement 
between the two simulations is excellent for all halo masses 
(excluding bins containing halos with fewer than 100 par- 
ticles). Combining the two allows for a consistent measure- 
ment of the halo mass function over seven decades in halo 
mass. At z = 6.2, the MS-II mass function lies systemati- 
cally above that of the MS. The most likely explanation of 
this difference is cosmic variance: the halos probed by ei- 
ther simulation at z = 6.2 are inherently rare objects, as the 
characteristic mass M* is 4.5 x 10 5 ft -1 Mp, at that tima^ 
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Figure 10. Halo bias at redshift zero. We combine results from 
the MS-II (filled circles) and the MS (open squares) to explore 
bias from 10 -4 to 10 M*. As expected, the bias decreases as the 
halo mass decreases, reaching b(M*) rs 1. At very low masses 
[Mv/M* <2x 10~ 2 or v < 0.55), the bias reaches an asymptotic 
value of 0.65. 



Furthermore, the MS-II probes only l/125 th the volume of 
the MS, making statistical fluctuations much more likely. 

In order to estimate the effects of cosmic variance on 
these mass functions, we divided the MS into 125 disjoint 
sub-cubes, each with the same volume as the MS-II, and 
we measured the scatter in mass functions and in the mean 
matter densities p m computed from these sub-volumes at 
z = 6.2. The full range of these mass functions is plotted as 
a gray shaded region in Figure[9] while the rms values at each 
mass are shown as black error bars on the MS data points. 
The MS-II points typically lie slightly outside of the rms 
region but well within the full distribution of mass functions, 
indicating that they are fully consistent with the MS when 
the volume of the MS-II is taken into account. We emphasize 
that the variation in the mass functions between the 125 
MS sub-cubes is not due to differences in the mean matter 
density, as the rms scatter in p m is only 2% while the rms 
scatter in the mass function exceeds 8% (the full range of 
the scatter exceeds ±20%) for all of the data points. 

4.2 Bias 

Dark matter halos do not cluster in the same way as the 
underlying mass density field but rather exhibit a bias rela- 



tive to the dark matter. Mo & White ( 1996 1, building on the 



earlier work of Efstathiou et al. ( 1988 1 and Cole & Kaiser 
( |1989[ ), showed that the two-point correlation function of 
halos should be simply related to that of the mass density 



10 The minimum halo mass in the MS, 10 10 h 1 Mq , corresponds 



to a peak height v = 8 c /cr(M, z) of 1.5 at z = 
equivalent to a mass of 7 X 10 13 h- 1 Mq at z = 0. 



6.2, which is 
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field. According to the excursion set model, on large scales 
one should find that 

ihh{M, z; r) = b 2 (M, z) £ mm (z; r) , (9) 

where the bias factor b is given by 

b(M,z) = l + S- 1 [y 2 -1). (10) 

Massive halos (M > M*) are therefore predicted to cluster 
more strongly than the underlying mass density field while 
low-mass halos should cluster less strongly. This basic pic- 
ture has been extensively validated, with newer models mak- 



ing improved quantitative predictions for bias (Jing 



1998 



Sheth fc Tormenl[T999l |Sheth et al.|[200l] |Seljak fc Warren| 



2004) 



The bias of low-mass halos (M <C M*) remains an un- 
resolved issue. In the Mo & White model, b -> 1 - S' 1 « 0.4 
for M <C M*.|Sheth fc Tormen| (|1999|) and|Seljak fc Warren 



(20041, on the other hand, find b ~ 0.7 in this same regime. 



The mass resolution of the MS-II allows us to study bias 
for M <C M*: halos with 200 particles, whose number den- 
sities and spatial distributions are certainly well-resolved, 
correspond to 2 x 10 -4 M t at z = 0. 

Halo bias at redshift zero is shown in Figure [10] for 
masses down to 2 x 10~ 4 A/*. The bias for M v < M + is 
constant at b ~ 0.65 over approximately two decades in 
mass, from M v /M* = 2 x 10~ 4 to 2 x 10" 2 . In terms of peak 



height v (shown on the upper horizontal axis of Figure 10 1, 
the bias is constant for v < 0.55. Above 0.1 M*, the bias 
rises rapidly with mass, reaching b = 1 at M slightly greater 
than M* and 6 » 1.5 for 10 M*. 

Figure[l0]also shows the predictions for b(M) from three 
fitting formulae: the original Mo & White prediction (solid 
curve), the Sheth & Tormen model (dashed), and the fit 
from Seljak & Warren (their eq. 5; dotted). The Seljak & 
Warren fit clearly agrees best with the data over the range 
plotted, though it slightly underpredicts the bias at M > M* 
and slightly overpredicts it at M < A/*. These differences 
are only at the 5% level, however. 



5 HALO FORMATION 
5.1 Formation Times 

The hierarchical nature of ACDM models means that the 
typical formation redshift Zf of halos with mass M is a de- 
creasing function of M. The form of the relation between Zf 
and M and its intrinsic scatter are important for a number 
of applications, such as understanding how well galaxy prop- 
erties can be predicted by halo mass alone, and how well a 
halo's history can be predicted from its present-day proper- 
ties. Such characterizations are complicated by the fact that 
the most useful definition of formation time depends on the 
question one is askinj 1 " 1 "] For example, the innermost region 
of a halo - where the galaxy resides - typically assembles 
much earlier than the outer regions which contain most of 
the mass ( |Zhao et al.|2003b| [Gao et al.|2004a l. 

One of the simplest definitions of formation redshift is 
the time at which a halo's main progenitor reaches a fixed 
fraction of its present-day mass. We use the merger trees 
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Figure 11. Mean half-mass formation redshift as a function of 
mass for halos from the MS-II (filled blue squares) and the MS 
(open red squares). We also show the relation for the 16% earliest 
and latest forming halos (dashed lines) and the best-fitting linear 
relation between log(l + Zf) and log Mi, (black dotted line; see 
Equation (JTTJl for the fitting parameters). This fit deviates from 
the mean relation by less than 2% over the entire range of masses 
plotted. 



described in Section 12.31 to trace each FOF halo back in 
time and define its formation redshift Zf as the first redshift 
at which one of the halo's progenitors reached half of the 
halo's redshift zero mass (we interpolate between snapshot 
redshifts to obtain Zf). This 'half-mass' formation time is the 
most common choice of formation time in the literature. We 
use M v as the definition of halo and progenitor mass when 
estimating such formation times (we have checked that the 
following results are insensitive to halo mass definition). 

In Figure [TT| we show the mean relation between halo 
mass M v and formation time Zf. In order to determine 
what mass is required for converged results, we compute 
Zf from the MS and compare with the MS-II. We find that 
the two simulations are in excellent agreement above a red- 
shift zero mass of 10 11 ft -1 Mq, corresponding to approx- 
imately 150 particles in the MS. This is the convergence 
limit we adopt, so we consider all halos with masses greater 
than 10 9 h' 1 Mq in the MS-II. We only include halos that 
SUBFIND determines to be bound, although in practice, this 
restriction makes almost no difference as the fraction of 
subhalos with N p > 100 that are unbound is very small. 
Over the entire range where halo formation can be resolved 
(9 < log 10 (M„/ft~ 1 M Q ) < 14.7), a simple linear fit in 
log(l + z/) versus log(M„) provides an excellent description 
of the data: 



1 + z f = 2, 



My 

10 10 h- 1 Mq 



(11) 



The maximum deviation between the binned data and Equa- 
tion (111 is 1.8% over the entire region where: (1) there are 
See, e.g., Li et al.|2008 for several possible definitions of Zf at least 100 halos per bin; and (2) halos have at least 125 
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Figure 12. Assembly bias from the MS-II (open symbols) and 
the MS (filled symbols). We plot the bias ol the oldest 20% of 
halos (red points) and the youngest 20% of halos (blue points), as 
well as the bias of the full halo sample (black points), as a function 
of halo mass. The oldest halos cluster much more strongly than 
the youngest, a ratio of 2.85 in the bias (which is over a factor 
of 8 in correlation amplitude). At low masses, the bias for each 
subset reaches an asymptotic value (si 0.4 for the young halos, 
Ri 1.15 for the old halos). 



particles. The 1-er scatter in the relation, defined as the log- 
arithmic difference between the 84 th and 16 th percentiles of 



the data (which are shown in dashed lines in Figure 111, de- 
creases gradually from <ri og i+ z = 0.6 to 0.4 as M v increases 
from 10 9 to 5 x 10 14 h' 1 M m . 



Neto et al. (20071 previously studied the relation be- 



tween median Zf and halo mass for a set of relaxed halos 
from the MS. Their fit is similar to ours although the pa- 
rameters differ slightly (their exponent is -0.046 and their 
normalization is approximately 2.74 using our form of the 
fitting formula) due presumably to the selection criteria used 
for the halos they studied, their use of M200 rather than 
M v , and the difference between the median and the mean 
relation. [McBride et al.| ( |2009[ ) have also recently computed 
the formation times of massive halos from the MS and fit- 
ted to the relation z f = alog 10 (M/10 12 M Q ) + b. This gives 
very similar results to ours over the range of the data they 
used (M > 10 12 h' 1 M©), with differences at the 5% level 
after an empirical normalization correction due to a differ- 
ent mass definition. Their formula underestimates 2/ from 



the MS-II at lower masses: 
approximately 20%. 



at 10 h Mq, the difference is 



5.2 Clustering and Formation Times 

In the simple version of the excursion set model of struc- 
ture formation, which is based on top-hat fc-space filtering, 
halo formation is governed by Markov random walks of the 
(linearly-extrapolated) mass density field <5(x;M). A direct 
consequence of this aspect of excursion set theory is the pre- 



diction that properties such as the clustering of halos depend 
on halo mass alone and not on halo assembly history. Recent 
Af-body simulations have produced results that contradict 
this prediction, however. Gao et al. (20051 used the MS to 



show that clustering depends strongly on formation time at 
masses M < M + : they found that early- forming halos clus- 
ter much more strongly than late-forming halos, indicating 
an "assembly bias". Subsequent work confirmed these find- 
ings and extended the results to the M > M* regime and to 
halo properties other than formation time, including concen- 
tration, substructure content, and spin (e.g., |Harker et al 
2006] [Wechsler et al.||2006] | Wetzel et al.||2007[ |Jing et aT 
2007||Bett et al.|2007|[Gao fc White|2007[ [Dalai et al.|2008] 



Angulo et al.|2008 l 

With the MS-II, we are able to investigate assembly 
bias at much lower masses than was previously possible: 
M ~ 10 -4 Mi, or equivalently v « 0.32. We split each mass 
bin into the oldest and youngest 20% of halos and com- 
pute the bias factor b(M) in the same manner described in 
Section [4. 2| Our results for the dependence of clustering on 
formation time are presented in Figure [12] We plot the bias 
of halos as a function of mass at redshift zero for the entire 
sample of halos (black symbols) as well as for the oldest 20% 
(red symbols) and youngest 20% (blue symbols) of halos at 
each mass. We show results for the MS-II (filled circles) for 
M < 2 x 10 11 h' 1 Mq and for the MS (open squares) at 
higher masses in order to maximize statistical significance. 

Over virtually the full range of masses probed by the 
MS-II, the biases of the oldest and youngest halos are ap- 
proximately constant. The bias of young halos is substan- 
tially lower than that of old halos, however, in agreement 
with previous work. The oldest 20% of halos at a given mass 
have a slight positive bias with respect to the dark matter 
distribution. The youngest 20% have a bias that is approx- 
imately b — 0.4; interestingly, and perhaps coincidentally, 
this is the value of 1 — S^ 1 ~ 0.4 predicted in the v < 1 



regime by Press-Schechter and excursion set models ( Cole & 
Mo fc White|1996| see also |Dalal et al 



Kaiser 1989 



Li et al. 



2008) 



(2008) have suggested that the magnitude of 
assembly bias depends on the adopted definition of for- 
mation time. To investigate whether our definition of for- 



mation redshift, M(z/) 



M(z = 0), influences our re- 



sults, we have repeated our analysis with a new definition: 
M(zf) = \ M(z = 0). We have also tested whether comput- 
ing the half-mass formation time relative to Mfof rather 
than to M v affects our results. Neither of these changes has 
any detectable influence, so the results we obtain using the 
standard half-mass formation time appear robust. 



6 'MILKY WAY' HALOS IN MILLENNIUM-II 

The Milky Way can provide us with a unique variety of in- 
sights into galaxy formation, so it is extremely interesting 
to study the formation of Milky Way-mass halos for com- 
parison with available and forthcoming observations of the 
detailed structure of our Galaxy. Cosmological simulations 
of representative volumes of the universe can provide large, 
statistically complete samples of Milky Way-mass halos, but 
they cannot resolve the full range of observable structures 
within the Milky Way. Even at the mass resolution of the 
MS-II, for example, the halos of dwarf spheroidal galaxies 
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167.18 
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Table 2. Comparison of Aquarius level 2 halos (Aq-X) and 
Aquarius halos in the MS-II (MS-II-X). M200771 is the mass of 
the dark matter halo, V max is the maximum value of the circular 
velocity curve, i? max is the radius at which the circular velocity 
curve attains its maximum, and V max su b is the circular velocity 
curve maximum for the largest subhalo in each halo; all values 
are quoted at 2 = 0. 



are just barely resolvable. An alternative tack is to focus all 
of one's computational resources on the formation of a single 
galaxy-mass halo, thus allowing substantially enhanced mass 
resolution piemand et al.||2007| |2008} |Springel et al.||2008 



|Stadel e"t~a l. 2008). With this method, one sacrifices statis- 
tical understanding of a representative sample for detailed 
analysis of one or a few objects. Here we discuss how these 
two approaches may be combined to extract maximum in- 
formation about the formation and evolution of galaxy-scale 
dark matter halos. 



6.1 The Aquarius Project and Millennium-II 

The Aquarius project ( Springel et al.|2008 I is a suite of cos- 
mological simulations of the formation of Milky Way-mass 
dark matter halos. Six halos (denoted 'Aq-A' through 'Aq- 
F') were simulated at up to five different levels of mass res- 
olution (levels 5 through 1). The highest resolution simu- 
lation, Aq-A-1, used a particle mass of 1.25 x 1O 3 /i _1 M , 
resulting in approximately 1.5 billion particles within Jfeoom 
at z = 0. The halos resimulated in the Aquarius project were 
select ed from the cosmological simulation 'hMS' (Gao et al. 



2008[ ), which followed 900* particles in a 100 h~ L Mpc box. 
Both the cosmological parameters of hMS and the ampli- 
tudes and phases of modes used to generate its initial con- 
ditions are identical to those in the MS-II. As a result, all 
the structures present in hMS are also present in the MS- 
II, but with a mass resolution that is a factor of 13.8 times 
better. Since the Aquarius halos were selected from hMS, 
they are also present in the MS-II and we can compare their 
properties in the MS-II and in the much higher resolution 
Aquarius resimulations. 
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Figure 13. Mass accretion histories for the six Aquarius level 
2 resolution halos (dashed lines) and for the corresponding MS- 
II halos (dotted lines). Even though the mass resolution of the 
Aquarius simulations is one thousand times better than that of 
the MS-II, the individual mass accretion histories of the Aquarius 
halos are captured quite well in the MS-II. 



6.2 Comparing Aquarius and Millennium-II halos 

As a first test, we compare some basic properties of the 
Aquarius halos in the Aquarius level 2 simulations - where 
the particle mass is (0.5 - 1) x 10 4 h' 1 M Q - and in the MS- 
II, where the mass resolution is approximately 1000 times 
lower. Table [2] contains such a comparison, with halos from 
the Aquarius simulations in the upper of each set of two rows 
(labeled 'Aq-X') and halos from the MS-II in the lower of 
each set of rows (labeled 'MS-II-X'). The first data column 
compares Afeoom values for the halos, showing that their 
masses agree very well: each M200™, measured from the MS- 
II agrees with the corresponding Aquarius resimulation to 
better than 6% and, for 3 of the 6, to better than 3%. The 
measured V ma .x values all agree within 3% and, for 5 of the 
6, to better than 1.5%. The radius at which the circular ve- 
locity curve peaks, -Rmax, typically agrees within 5%, with a 
maximum deviation of 10%. Table [5] also lists the maximum 
circular velocity of the largest subhalo in each simulation. 
In all cases this is the same subhalo in the MS-II and the 
Aquarius resimulation, and the circular velocities typically 
agree to within a few percent. 

A more stringent test is to consider the mass accretion 
histor^p^for each Aquarius halo and to compare results from 
the Aquarius resimulations and from the MS-II. For the MS- 
II Aquarius halos, we use the merger trees described in Sec- 
tion |2.3| to determine the main progenitor of each halo at 
each redshift, and we define the mass accretion history of a 



12 for a discussion of the statistical properties of mass accretion 
histories, see |Lacey fc C ole (1993); Wechsler et al. (2002); v an den| 
^( p?002j l; |Zhao et al.| 



Bosch 



<|2003bja| [20081: 



jMcBride et al.| ( |2009| . 
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halo as the mass of its main progenitor M(z) at each red- 
shift z. Merger trees have also been built for the halos in 
the Aquarius resimulations, and we use these in an identical 
manner to define the corresponding mass accretion histories. 

The results are compared in Figure |13| where we plot 
mass accretion histories over the redshift range < z < 9. 
We use dashed lines for the Aquarius resimulations and dot- 
ted lines for the MSTI. The assembly histories of all the 
halos are reproduced remarkably well at MS-II resolution. 
This is a highly non-trivial test, as the mass resolution of 
the Aquarius level 2 simulations is one thousand times bet- 
ter than that of the MS-II: at z — the Aquarius level 2 
halos have on the order of 1.80 x 10 8 particles within 7?200m 
while their MS-II counterparts have approximately 1.80x 10 5 
particles within this radius. 

This agreement is a testament to the integration accu- 
racy of the GADGET-3 code and shows that the properties 
of Milky Way-mass halos and their most massive subhalos 
are well converged in the MS-II. As a result, the MS-II will 
be very useful for understanding the statistical properties of 
Milky Way-mass halos since it contains over six thousand 
halos with 11.5 < log 10 (M„//i _1 M Q ) < 12.5. These halos 
can be used for a detailed statistical investigation of the 
growth, internal structure and subhalo populations of halos 
similar in mass to that of our own Galaxy (Boylan-Kolchin 
2009, in preparation). This will test the extent to which the 
six Aquarius halos are representative of the full halo popu- 
lation at this mass, thereby allowing results obtained from 
the ultra-high resolution Aquarius resimulations to be inter- 
preted in their full cosmological context. 



7 DISCUSSION 

Understanding galaxy formation in a ACDM universe re- 
quires knowledge of physical processes over a wide range of 
scales, from sub-galactic to cosmological. Simulations with 
volumes large enough to probe the statistics of large-scale 
structure and resolution high enough to resolve subhalo dy- 
namics within galaxy halos are thus critical for this quest. 
We have presented initial results from a new simulation, the 
Millennium-II Simulation, that can resolve all halos down to 
mass scales comparable of the halos of Local Group dwarf 
spheroidal galaxies. 

Furthermore, the Millennium-II Simulation is closely 
connected to two other very large computational endeav- 
ors, the Millennium Simulation and the Aquarius Project. 
Throughout this paper we have shown that MS-II results 
agree extremely well with those from the MS, so we can 
combine the two simulations to cover an even broader range 
of physical scales. We have also shown that the properties of 
the resimulated Aquarius halos agree precisely and in con- 
siderable detail with those of their counterparts in the MS-II, 
even though the resimulations have 1000 times better mass 
resolution. This gives us confidence not only that the assem- 
bly histories of Milky Way-mass halos are well resolved in 
the MS-II, but also that the properties of their more massive 
subhalos are converged. As a result, we will be able to use 
the MS-II to make statistical statements about an ensemble 
of galaxy-scale halos where the Aquarius Project provides 
much higher resolution results for a limited but well under- 
stood subset. 



There is much to do with the MS-II, both from a dark 
matter perspective and from the point of view of galaxy 
formation. There are four cluster-size halos in the MS-II 
with over 60 million particles at z — 0; investigating den- 
sity profiles and substructure abundances for these objects, 
and comparing with state-of-the-art galaxy-scale simula- 
tions such as Aquarius, will shed light on whether dark mat- 
ter structures are self-similar as a function of scale. Subhalo 
survival times and merger rates, which are crucial ingredi- 
ents in galaxy formation models, can be checked in untested 
regimes. Furthermore, directly resolving much smaller halos 
means that semi-analytic models initially developed for the 



Millennium Simulation 


Springel et al. 2005| Croton et al. 


2006 Bower et al.|2006 


De Lucia & Blaizot|2007 1 can now 



be updated and extended to much lower galaxy masses (Guo 
et al. 2009, in preparation). With this paper, we publicly re- 
lease the FOF halo and subhalo catalogs and the subhalo 
merger trees in a searchable databasq_j structured in the 
same way as has already been done for the Millennium Sim- 
ulation (Lemson & The Virgo Consortium 2006). As work 



progresses, we plan to make MS-II semi-analytic galaxy cat- 
alogs available as well. This will allow others to use the MS-II 
data for their own research purposes. 
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